Tutorial Part IV — Getting started with animal movement

Introduction

In this tutorial, I will show you how to get started with (real) animal movement data. To put it simple, animal movement data are just points in time and space, i.e. re-locations of an individual that can be provided as a data.frame with four columns. An identifier for the individual is only needed if you have various animals in one data set, i.e. if you have collared several individuals with the same collar or already appended several data sets:

Identifier Timestamp x_coordinate y_coordinate
Fox_1 2022-01-01 08:10:04 13.47027 52.497802
Fox_1 2022-01-01 08:15:04 13.47015 52.497813

From Course2_SpatialR, you might remember that such data sets can be imported as simple .txt or .csv files. Once you know in which coordinate reference system (CRS) your coordinates are stemming from, you can assign it to the data.frame, thereby you are creating a geo-referenced data set: a SpatialPointsDataframe-Object (if you work with the R-package {sp}) or an sf-Object (with package {sf}). Looking at these x- and y-coordinates, you might remember that they look like angular units, e.g. decimal degrees, and the x-coordinate refers to longitude and the y-coordinate to latitude. Hence, we are dealing with EPSG-code 4326. Once you have imported and geo-referenced your data, you can plot and thoroughly check the data.

Checking your data is a mandatory prerequisite before any analysis (see Zuur et al. 2010). With your movement data, you might have outliers in space (e.g. when you test your collars in location A which might be hundreds of kilometers from your actual trapping site), your data might not have been regularly sampled because the animal was hiding and there was no access to the satellite (weak GPS-signal) etc etc. This uneven sampling can affect calculations of speed, for example, and therefore a series of packages have been developed to deal with those issues in movement data, like {move}.

The most prominent movement analyses comprise home range estimation, calculations of habitat preference, behavioral analyses (e.g. Hertel et al. 2021) and detection of movement syndromes (personality differences) (Michelangeli et al. 2021). Some of the packages listed in the next chapter are designed for these analyses.

In the following, we will step by step start with loading and exploring movement data of the female red fox (Vulpes vulpes) ‘Q von Stralau’, who had a small home range in Berlin, Germany. She was collared in January 2018 and had a very stable daily routine as shown by 4 (20) min relocation intervals. Her location data are stored in the file with the tag number tag5334_gps.txt. For any details on the data please refer to Kimmig (2021) and Scholz (2020).

We will, however, only use the data of one month in this exercise, as data sets quickly get too big. Please note: the courses Course1_IntroR and Course2_SpatialR are obligatory for this tutorial.

Useful (web)sites and reading

  • For analysis of telemetry data: packages adehabitatHR, adehabitatLT, move,recurse, momentuHMM, moveHMM, ctmm, amt

Methods papers

Example papers

Getting started

To follow the tutorial, you can either clone or download the repository at https://github.com/stephkramer/Course4_MoveQ or you create your own R-project, copy the raw data and type the code chunks into an R-Script. Please refer to the section on using R-projects in Course2_RSpatial.

If you start with your own R-project, I strongly recommend to use the {d6}-package Cedric Scherer provided. This package automatically sets up the ideal folder structure: https://github.com/EcoDynIZW/d6

In any case, my course folder has the following structure:

.
└── Course4_MoveQ         – (root folder)
    ├─── data-raw         – (contains the file with the GPS locations)
    ├─── docs
    ├─── output           – (contains the cleaned files and processed data)
    ├─── plots
    ├─── R                – (contains the R-script)
    └ Course4_MoveQ.Rproj – (the RStudio-Project file location)


Necessary packages to install and load

We first have to install the packages and load them before we can use the functions we need.

## package names:
# pkgs = c("here", "sp", "sf", "dismo", "raster", "GISTools",  "rgdal", 
#          "maptools", "rgeos", "rgl", "rasterVis", "adehabitatHR", "move", "tmap",
#          "plotly", "circular", "gganimate", "moveVis", "ggmap", "maps", "mapproj",
#          "viridis", "dplyr", "devtools", "lubridate", "patchwork", "gower",
#          "colorspace", "ragg", "ggtext", "pdftools", "units", "leaflet",
#          "glue", "cowplot", "tidyverse", "ggplot2", "Cairo") 
pkgs = c("here", "lubridate", "sf", "sp", "raster", "rgeos", "ggplot2", "tmap",
         "circular", "move", "adehabitatHR", "viridis", "devtools", "plotly") 

# install.packages(pkgs) # only run this line once for installing the packages!
# update.packages()

Tip of the day: If you have already installed some of the packages above, you can first check which ones are already installed, and save the ones not installed in an object called ‘my_packages’ and only install the missing ones:

my_packages <- pkgs[!(pkgs %in% installed.packages()[,"Package"])]
my_packages
## character(0)
if(length(my_packages) > 0) install.packages(my_packages)

Now load the packages:

library(here)          ## for easy directory management
library(lubridate)     ## for date handling
library(sf)            ## for handling simple feature objects
library(sp)            ## for handling Spatial* objects
library(raster)        ## for handling raster objects
library(rgeos)         ## retiring end of 2023 -> use sf
library(ggplot2)       ## for nice static plots
library(tmap)          ## for nice interactive maps
library(circular)      ## for circular stats + plots
library(move)          ## for `move` objects
library(adehabitatHR)  ## adehabitatLT  #adehabitatMA                           
library(viridis)       ## for perceptually uniform color palettes
# library(dismo)
# library(GISTools)
# library(rgdal)    # retiring end of 2023 -> use stars/ terra /sf
# library(maptools) # retiring end of 2023 -> use sf
# library(rgl)
# library(rasterVis)
# library(plotly)
# library(dplyr)
# library(gganimate)
# library(moveVis)
# library(ggmap)
# library(maps)  
# library(mapproj)
# library(devtools)
# library(glue)
# library(cowplot)
# library(tidyverse)
# library(Cairo)
# library(colorspace)
# library(ragg)
# library(ggtext)
# library(pdftools)
# library(units)
# library(patchwork) # to combine plots
# library(leaflet)
# library(gower)

Set the working environment

Now that we are working inside an R-project, we can use the easy functionality of the {here} package, which is automatically pointing to your project’s root folder, e.g.:

here::here()

Hence, there is no need to use the function setwd() any more. Note: if it does not work, please close RStudio, go to your Explorer and double-click on the .Rproj file. Then, under ‘files’ (usually lower right panel) double-click on the R folder and open the script.

Load data

The movement data are stored in the data-raw subfolder. Let’s check which files are available.

lf <- list.files(path = here("data-raw"), full.names = TRUE) 
lf
## [1] "C:/Users/admin/Course4_MoveQ/data-raw/DieTiere_32633.cpg"
## [2] "C:/Users/admin/Course4_MoveQ/data-raw/DieTiere_32633.dbf"
## [3] "C:/Users/admin/Course4_MoveQ/data-raw/DieTiere_32633.prj"
## [4] "C:/Users/admin/Course4_MoveQ/data-raw/DieTiere_32633.qpj"
## [5] "C:/Users/admin/Course4_MoveQ/data-raw/DieTiere_32633.shp"
## [6] "C:/Users/admin/Course4_MoveQ/data-raw/DieTiere_32633.shx"
## [7] "C:/Users/admin/Course4_MoveQ/data-raw/geo-raw"           
## [8] "C:/Users/admin/Course4_MoveQ/data-raw/KettleHoles.txt"   
## [9] "C:/Users/admin/Course4_MoveQ/data-raw/tag5334_gps.txt"

The output lists two results, there is another subfolder, and many files. E.g., the element 7 of the vector lf, lf[7], is a folder. If you already know that your movement data contains e.g. ‘gps’ in its name or is stored as ‘.txt’ or ‘.csv’ files, you can directly search for those files with the pattern argument:

## check the difference, and note: full.names is set to FALSE
thefile <- list.files(path = here("data-raw"), pattern = "gps", full.names = FALSE)
thefile
## [1] "tag5334_gps.txt"
## ...and here to TRUE
thefullfile <- list.files(path = here("data-raw"), pattern = "gps", full.names = TRUE)
thefullfile
## [1] "C:/Users/admin/Course4_MoveQ/data-raw/tag5334_gps.txt"

Now load the data file. Our first fox filename should be tag5334_gps.txt.

## Let's get the number of the list in lf that corresponds with the filename.
## Tip of the day: always code everything you can, do not insert numbers or filenames
## by hand, i.e. lf[2]. If you for example add new data, the numbering changes!
lf_number <- which(lf == thefullfile)
dat_anim <- read.table(file = lf[lf_number], header = TRUE, fill = TRUE, sep = ",")

## or, since with 'thefullfile' you gave the complete path, you can also directly load it
# dat_anim <- read.table(file=thefullfile, header=TRUE, fill=TRUE, sep=',') ## CED: better use this, the other one is pretty confusing for them? Or do you want to show it so they can access multiple files with this approach later?


Data check and cleaning

Checking for missing or incomplete information

Let’s have a look at the data. This is the typical way data are stored on e-obs collars:

dat_anim[1:5,] ## recap: head(dat_anim) also works
##   key.bin.checksum tag.serial.number         start.timestamp longitude latitude
## 1       2135768252              5334 2018-10-30 04:00:00.000  13.47047 52.49499
## 2       2210577241              5334 2018-10-30 04:20:00.000  13.47157 52.49526
## 3       1824790132              5334 2018-10-30 04:40:00.000  13.47103 52.49525
## 4       4079427513              5334 2018-10-30 07:00:00.000  13.47109 52.49526
## 5        380027806              5334 2018-10-30 08:00:00.000  13.47105 52.49527
##   height.above.ellipsoid type.of.fix status used.time.to.get.fix
## 1                   76.3           3      A                   22
## 2                   81.6           3      A                   18
## 3                   66.5           3      A                   53
## 4                   73.5           3      A                   24
## 5                   71.1           3      A                   47
##          timestamp.of.fix battery.voltage fix.battery.voltage temperature
## 1 2018-10-30 04:00:23.000            3604                3389          10
## 2 2018-10-30 04:20:19.000            3616                3333          11
## 3 2018-10-30 04:40:54.000            3624                3374          13
## 4 2018-10-30 07:00:25.000            3659                3491          11
## 5 2018-10-30 08:00:48.000            3660                3470          13
##   speed.over.ground heading.degree speed.accuracy.estimate
## 1              1.64         233.31                    1.67
## 2              4.43         327.02                    3.79
## 3              0.07           0.00                    0.73
## 4              0.01           0.00                    0.65
## 5              0.04           0.00                    0.89
##   horizontal.accuracy.estimate
## 1                         7.94
## 2                        59.65
## 3                        15.36
## 4                         7.17
## 5                         8.45

For now, we will only work with few columns. tag.serial.number refers to the individual identifier (collar ID). There are two columns with timestamps: start.timestamp is the preprogrammed time-interval. Then there is the timestamp.of.fix, which is the real time the GPS-location was recorded. This is usually a bit later, i.e.
( = start.timestamp + used.time.to.get.fix + 1 second), as it takes some time for the collar unit to connect to the satellite. The spatial info is stored in the columns longitude and latitude.

Before you can transform the data.frame dat_anim into a georeferenced spatial object, you need to check whether there are missing locations in your data.frame, otherwise you will get an error message on transformation.

This can happen if no GPS-signal could be recorded. Or - importantly - some collars are only activated when the animal is moving to save battery life. In that case, the missing GPS coordinates would correspond with the last position (= be the same). Depending on which analysis you want to do, e.g. define resting places, you might need to fill the missing positions again.

In our case, there are a lot of missing values in the locations:

## if the latitude-entry is missing, the longitude value will also be missing
## so it is enough to only check the latitude
which(is.na(dat_anim$latitude))
##   [1]   48   65   89   90   91  129  183  185  222  224  226  229  260  266  355
##  [16]  358  415  421  422  423  457  474  510  511  577  601  616  621  626  646
##  [31]  673  674  685  696  721  774  784  796  822  833  842  855  863  907  921
##  [46]  940  993 1006 1038 1070 1077 1078 1088 1126 1150 1159 1183 1188 1215 1238
##  [61] 1268 1272 1304 1330 1480 1516 1550 1552 1559 1600 1637 1680 1717 1721 1736
##  [76] 1738 1742 1772 1820 1824 1830 1867 1884 1973 1997 1999 2020 2023 2024 2025
##  [91] 2026 2064 2075 2076 2109 2162 2234 2313 2329 2374 2379 2396 2419 2453 2455
## [106] 2518 2551 2609 2614 2638 2687 2728 2745 2770 2776 2783 2800 2804 2805 2854
## [121] 2926 2927 2931 2935 2943 2950 2961 2963 2968 2990 2991 3043 3044 3052 3090
## [136] 3102 3114 3172 3225 3226 3305 3328 3361 3372 3373 3386 3409 3498 3516 3517
## [151] 3518 3520 3591 3622 3623 3624

Delete rows with missing spatial info:

dat_anim_na <- dat_anim[!is.na(dat_anim$latitude),] ## alternatively, use complete.cases()

Make the crosscheck if there is missing info in a row in longitude. This should NOT be the case after we had deleted those rows:

which(is.na(dat_anim_na$longitude)) # none
## integer(0)

Checking for coarse spatial outliers

It could happen the collar was tested e.g. in Berlin, but the animal was finally caught and collared far away. Make a quick check whether there are strange locations:

plot(dat_anim_na$latitude ~ dat_anim_na$longitude)

There do not seem to be coarse outliers, data look compact.

Appending sunrise and sunset as columns (working with dates)

We will now add some additional columns, where separate days (numbered from 1 to 365), the month (from 1 to 12) and the hour of the day are stored (from 0 to 23). Sunset and sunrise can be calculated based on dates and the location (latitude, longitude). This can be done with the package {lubridate}. Check the vignette!

## have a look at the timestamp format, here the first row of the column start.timestamp
dat_anim_na$start.timestamp[1]
## [1] "2018-10-30 04:00:00.000"
## define date-time format - the format is year-month-day_hour:min:sec:
dat_anim_na$start.timestamp <- ymd_hms(dat_anim_na$start.timestamp, tz="Europe/Berlin") 
dat_anim_na$start.timestamp[1]
## [1] "2018-10-30 04:00:00 CET"

Now append the information:

dat_anim_na$yearday <- yday(dat_anim_na$start.timestamp)
dat_anim_na$month   <- month(dat_anim_na$start.timestamp)
dat_anim_na$hour    <- hour(dat_anim_na$start.timestamp)
dat_anim_na$kweek   <- week(dat_anim_na$start.timestamp)
dat_anim_na$date    <- date(dat_anim_na$start.timestamp)
## crosscheck with
# head(dat_anim_na)

In addition, we want to calculate the hours of sunset and sunrise as well as daylength. For this, we need to install a package that is still under development, i.e. which is not on CRAN. We therefore must download and install it locally:

# devtools::install_github("bgctw/solartime") ## run only once to install package from GitHub
library(solartime)

For computing sunset and sunrise, the latitude and longitude must be provided as well:

dat_anim_na$sunrise   <- computeSunriseHour(timestamp = dat_anim_na$start.timestamp,
                                            latDeg = dat_anim_na$latitude,
                                            longDeg = dat_anim_na$longitude)

dat_anim_na$sunset    <- computeSunsetHour(dat_anim_na$start.timestamp,
                                           dat_anim_na$latitude,
                                           dat_anim_na$longitude)

dat_anim_na$daylength <- computeDayLength(dat_anim_na$start.timestamp,dat_anim_na$latitude)

dat_anim_na$daytime   <- computeIsDayByLocation(dat_anim_na$start.timestamp,
                                                dat_anim_na$latitude,
                                                dat_anim_na$longitude)

## check new variables
# head(dat_anim_na)
hist(dat_anim_na$daylength)

unique(dat_anim_na$daytime)
## [1] FALSE  TRUE

Check for temporal outliers

Check if there are strange dates, or dates before you collared the animal (e.g., usually the collar is activated and tested before the animal is collared):

table(dat_anim_na$date)
## 
## 2018-10-30 2018-10-31 2018-11-01 2018-11-02 2018-11-03 2018-11-04 2018-11-05 
##         24         40         35         39         42         34         34 
## 2018-11-06 2018-11-07 2018-11-08 2018-11-09 2018-11-10 2018-11-11 2018-11-12 
##         41         42         42         41         38         38         39 
## 2018-11-13 2018-11-14 2018-11-15 2018-11-16 2018-11-17 2018-11-18 2018-11-19 
##         41         38         34         38         44         34         35 
## 2018-11-20 2018-11-21 2018-11-22 2018-11-23 2018-11-24 2018-11-25 2018-11-26 
##         36         34         40         46         40         37         41 
## 2018-11-27 2018-11-28 2018-11-29 2018-11-30 2018-12-01 2018-12-02 2018-12-03 
##         43         39         34         42         45         37         39 
## 2018-12-04 2018-12-05 2018-12-06 2018-12-07 2018-12-08 2018-12-09 2018-12-10 
##         38         39         38         31         38         36         36 
## 2018-12-11 2018-12-12 2018-12-13 2018-12-14 2018-12-15 2018-12-16 2018-12-17 
##         39         34         38         37         37         43         43 
## 2018-12-18 2018-12-19 2018-12-20 2018-12-21 2018-12-22 2018-12-23 2018-12-24 
##         37         31         41         36         38         42         44 
## 2018-12-25 2018-12-26 2018-12-27 2018-12-28 2018-12-29 2018-12-30 2018-12-31 
##         32         30         38         38         36         37         39 
## 2019-01-01 2019-01-02 2019-01-03 2019-01-04 2019-01-05 2019-01-06 2019-01-07 
##         33         43         37         36         39         39         45 
## 2019-01-08 2019-01-09 2019-01-10 2019-01-11 2019-01-12 2019-01-13 2019-01-14 
##         36         38         41         33         26         35         40 
## 2019-01-15 2019-01-16 2019-01-17 2019-01-18 2019-01-19 2019-01-20 2019-01-21 
##         28         35         35         40         39         34         31 
## 2019-01-22 2019-01-23 2019-01-24 2019-01-25 2019-01-26 2019-01-27 2019-01-28 
##         36         36         34         36         43         34         33 
## 2019-01-29 2019-01-30 2019-01-31 2025-12-26 
##         35         34         27          1

It might be easier to spot visually. We use the {ggplot2} package here as it recognises and repsects the date format:

ggplot(dat_anim_na, aes(date)) +
  geom_bar() +
  theme_bw()

## compare with plot(table(dat_anim_na$date))

There is a strange date → 2025-12-26!

Delete this data row:

delme <- which(dat_anim_na$date == "2025-12-26")
dat_anim_na[delme,]                 ## check observation
##      key.bin.checksum tag.serial.number     start.timestamp longitude latitude
## 2622       2948730959              5334 2025-12-26 04:00:00  13.47851 52.49133
##      height.above.ellipsoid type.of.fix status used.time.to.get.fix
## 2622                   78.8           3      A                   23
##             timestamp.of.fix battery.voltage fix.battery.voltage temperature
## 2622 2025-12-26 04:00:24.000            3582                3287           3
##      speed.over.ground heading.degree speed.accuracy.estimate
## 2622              0.03              0                    0.64
##      horizontal.accuracy.estimate yearday month hour kweek       date  sunrise
## 2622                         4.35     360    12    4    52 2025-12-26 8.390617
##        sunset daylength daytime
## 2622 15.80777  7.417151   FALSE
dat_anim_na <- dat_anim_na[-delme,] ## delete the strange date and 
table(dat_anim_na$date)             ## check again
## 
## 2018-10-30 2018-10-31 2018-11-01 2018-11-02 2018-11-03 2018-11-04 2018-11-05 
##         24         40         35         39         42         34         34 
## 2018-11-06 2018-11-07 2018-11-08 2018-11-09 2018-11-10 2018-11-11 2018-11-12 
##         41         42         42         41         38         38         39 
## 2018-11-13 2018-11-14 2018-11-15 2018-11-16 2018-11-17 2018-11-18 2018-11-19 
##         41         38         34         38         44         34         35 
## 2018-11-20 2018-11-21 2018-11-22 2018-11-23 2018-11-24 2018-11-25 2018-11-26 
##         36         34         40         46         40         37         41 
## 2018-11-27 2018-11-28 2018-11-29 2018-11-30 2018-12-01 2018-12-02 2018-12-03 
##         43         39         34         42         45         37         39 
## 2018-12-04 2018-12-05 2018-12-06 2018-12-07 2018-12-08 2018-12-09 2018-12-10 
##         38         39         38         31         38         36         36 
## 2018-12-11 2018-12-12 2018-12-13 2018-12-14 2018-12-15 2018-12-16 2018-12-17 
##         39         34         38         37         37         43         43 
## 2018-12-18 2018-12-19 2018-12-20 2018-12-21 2018-12-22 2018-12-23 2018-12-24 
##         37         31         41         36         38         42         44 
## 2018-12-25 2018-12-26 2018-12-27 2018-12-28 2018-12-29 2018-12-30 2018-12-31 
##         32         30         38         38         36         37         39 
## 2019-01-01 2019-01-02 2019-01-03 2019-01-04 2019-01-05 2019-01-06 2019-01-07 
##         33         43         37         36         39         39         45 
## 2019-01-08 2019-01-09 2019-01-10 2019-01-11 2019-01-12 2019-01-13 2019-01-14 
##         36         38         41         33         26         35         40 
## 2019-01-15 2019-01-16 2019-01-17 2019-01-18 2019-01-19 2019-01-20 2019-01-21 
##         28         35         35         40         39         34         31 
## 2019-01-22 2019-01-23 2019-01-24 2019-01-25 2019-01-26 2019-01-27 2019-01-28 
##         36         36         34         36         43         34         33 
## 2019-01-29 2019-01-30 2019-01-31 
##         35         34         27
plot(table(dat_anim_na$date))       ## plot the number of fixes per day

Save the cleaned file

Finally, we save the processed data file that we will use for exploration and analysis into the subfolder .../output/data-proc. There are two options we can use:

  • R data file - can only be opened/ read with R by using function readRDS()
saveRDS(dat_anim_na, file = here("output", "data-proc", "tag5334_gps_proc.Rds"))
  • interchange file format .csv
write.csv(dat_anim_na, file = here("output", "data-proc", "tag5334_gps_proc.csv"))

Check your output folder for these files. The efficient ‘.Rds’ file storage is about 3 times smaller than the ‘.csv’ file.

Exercise 1

This is a recap from Course2. Please plot the animal relocations in two different colours based on the column daytime, using the data.frame and one of the options to create maps with {ggplot2}, {ggpmap}, {leaflet} or {tmap}. Use the processed file of fox Q (tag5334_gps_proc) and do it in a separate script. Save your script as Course4_Exercise1_*yourname*.R.

Hint: Remember to load the relevant libraries. Note that you might want to plot the locations in the correct spatial dimensions by projecting it using the functions st_as_sf() and st_transform().



Data exploration

Load the cleaned data file

I recommend that you store the raw file safely and continue with the cleaned and processed file after major data manipulations have been conducted to minimize errors. I’d even suggest to make these steps in different R-scripts.

anim_proc <- readRDS(file = here("output", "data-proc", "tag5334_gps_proc.Rds"))
# head(anim_proc)

Have a look whether there are gaps/ missing days in the data, or whether we have approximately a regular number of fixes each day.

plot(table(anim_proc$date))

Again, we can use {ggplot2} to plot true dates:

ggplot(dat_anim_na, aes(date)) +
  geom_bar() +
  theme_bw()

Now let’splot the Julian day, yearday:

plot(table(anim_proc$yearday))  ## what is the difference?

The number of fixes (locations taken) per day seems to be quite regular. Mind the difference when plotting per date (sorted) and yearday (1-365).

Plot activity

Now we can have a look at the distribution of the fixes during the day, which is a hint on the active phase of the animal:

plot(table(anim_proc$hour))   ## plot the number of fixes per hour

More fixes during the night - this is the active phase. As these are circular data, we can plot the number of fixes as sign for activity (knowing that our collars did not record when foxes were inactive! Otherwise, we cannot distinguish missing data from inactive times! -> bias).

## simple circular plot (rose diagram)
timetoplot <- circular(anim_proc$hour %% 24, ## convert to 24 hrs = bins
                       units = "hours", template = "clock24")

## Note: use namespace `circular::` here as there are multiple functions called `rose.diag()`
circular::rose.diag(timetoplot, bin = 24, col = "blue",
                    main = "Events by Hour (sqrt scale)", prop = 3)

Check whether the activity during the day was before sunrise (our ‘daytime’ column which is either TRUE (=yes, daylight) or FALSE (dark hours):

# with ggplot
# code adapted from https://gist.github.com/mattbaggott/4361381
ggplot(anim_proc, aes(x = hour,fill = daytime)) +
  geom_bar(width = 1, color = "grey20") +
  coord_polar(start = -0.15) + ## to center setchange start
  theme_minimal() +
  scale_x_continuous(breaks = 0:23) +
  scale_fill_brewer(palette = "Set2", name = NULL, labels = c("Night", "Day")) +
  labs(x = NULL, y = "Count", title = "Events by Time of the Day")

Apart from the outliers between 11 and 13 o’clock, the active phase was during the dark hours.

Let’s fix the yearly split (we could also work with dates but take a manual route here):

Convert the data into a spatial object

Now we transform the data into sf and sp objects and assign the coordinate reference system:

# transform into spatial simple feature sf object
mydf_sf <- st_as_sf(x = data.frame(anim_proc),
                       coords = c("longitude", "latitude"),
                       crs = 4326,
                       sf_column_name = "geometry" )

# transform into SpatialPointsDataFrame  - for crosschecking
mydf_sp <- as(mydf_sf, "Spatial")

And then we project the reference system from angular units to a planar coordinate reference system in meters:

# transform CRS to projected one in meter distance units
mydf_sf_trans <-  st_transform(mydf_sf, 3035 )  # EPSG-code  
mydf_sp_trans <-  spTransform(mydf_sp, CRS("+init=epsg:3035"))

Recently, there are issues with the missing datum in the CRS-specifications. We ignore this for now. More on the issue of moving from proj4 to proj6 in the future:
* https://inbo.github.io/tutorials/tutorials/spatial_crs_coding/
* https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html

Make a quick plot of the data

Recap for styling the plot:
https://www.r-bloggers.com/2021/12/introduction-to-geospatial-visualization-with-the-tmap-package/

tmap_mode(mode = "view")

tm_shape(shp = mydf_sf_trans) + 
  tm_dots(size = 0.01, 
          col = "daytime",  
          alpha = 0.5)

Have a deep look at the data: did the fox really swim? (points in Lake Rummelsburg). Or are these outliers? No: the lake was frozen - check the date! Note: you really need to know your animals and the area before analysing data!

Basic movement metrics

Transform your data into a {move} object

check here for the possibilities of what to do with the package:
* https://cran.r-project.org/web/packages/move/vignettes/move.html
* https://rdrr.io/cran/move/f/vignettes/move.Rmd 

####### TODO use transformed animal data for spatial measures! #####
####### WTF funktioniert das nicht???????   ########
# Ich wollte andere Koordinaten nehmen, weil unten hrBootstrap 
# komische measures ausspuckt - viel zu kleine HRs

#mydf_sf_trans$coords <- st_coordinates(mydf_sf_trans)

## CED: Coords gibt einen Vekor aus. Der ist zwar als Coords.X und Coords.Y gelistet, allerdings nur über Coords[,1] und Coords[,2] ansprechbar. Also entweder einzeln zuweisen (wahrscheinlich leichter verständlich) oder im move() mit den eckigen Klammern referenzieren

## assign columns with coordinates to use later in data frame
mydf_sf_trans$long <- st_coordinates(mydf_sf_trans)[,1]
mydf_sf_trans$lat <- st_coordinates(mydf_sf_trans)[,2]

## new object as data frame
tt <- st_drop_geometry(mydf_sf_trans) ## CED: more descriptive name needed

## now turn into move object
my_fox <- move(x = tt$long, y = tt$lat,
               time = as.POSIXct(tt$start.timestamp, format = "%Y-%m-%d %H:%M:%S"),
               proj = CRS("+init=epsg:3035"),
               data = tt,
               animal = "FoxQvonStralau")  
## CED: nehme an, das kann dann weg?
# my_fox <- move(x = anim_proc$longitude, y = anim_proc$latitude,
#                time=as.POSIXct(anim_proc$start.timestamp,format="%Y-%m-%d %H:%M:%S"),
#                proj= CRS("+init=epsg:4326"),
#                data=anim_proc,
#                animal="FoxQvonStralau")  
# head(my_fox); nrow(my_fox)

Time lag, turning angle and speed

Time lags between the fixes:

timeLag(my_fox, unit = "mins")[1:5] ## varies a lot, but minimum here 20 min
## [1]  20  20 140  60 240
min(timeLag(my_fox, units = "mins"))
## [1] 17.25
mean(timeLag(my_fox, units = "mins"))
## [1] 38.48087
max(timeLag(my_fox, units = "mins"))
## [1] 240

Turning angles (makes only sense with high resolution data):

turnang <- angle(my_fox)
# turnang <- turnAngleGc(my_fox) # the same
# using the absolute abs(), because -180 and 180 has similar meaning:
# the animal keeps the direction
hist(abs(turnang))

Speed and step length:

steplength <- distance(my_fox) ## know your units!
hist(steplength)

max(steplength) ## max 1 km in 20 min = 3 km / h
## [1] 1028.777
hist(speed(my_fox)) ## units? -> check ?speed

The homerange concept - MCP, kernel, aKDE

See lecture!

Minimum Convex Polygon

Unfortunately, we have to use SpatialDataFrame-Objects for activity range calculation, as the package {adehabitatHR} is based on the old sp format.

You can calculate the total area covered, or, if you specify a column, you can calculate the activity range per daytime, animalID or months… Use the SpatialDataFrame-Object from the {sp} package we created above.

mcp_daytime <- mcp(mydf_sp_trans[,'daytime'], percent = 95, unout = "km2") # MCP (95%) 
mcp_daytime ## area for both (daytime separated) polygons
## Object of class "SpatialPolygonsDataFrame" (package sp):
## 
## Number of SpatialPolygons:  2
## 
## Variables measured:
##          id      area
## FALSE FALSE 0.5568531
## TRUE   TRUE 0.0599125
hrs <- mcp.area(mydf_sp_trans[,'daytime'], percent = seq(50, 100, by = 5), unout = "km2")

hrs ## home range size in km2 for the different MCP-levels
##         FALSE        TRUE
## 50  0.1600393 0.003554521
## 55  0.1748087 0.003672804
## 60  0.1974765 0.003801935
## 65  0.2233283 0.003978713
## 70  0.2573820 0.004021880
## 75  0.2860970 0.004097756
## 80  0.3632446 0.004344078
## 85  0.4269748 0.004867771
## 90  0.4988853 0.010032773
## 95  0.5568531 0.059912502
## 100 0.6841272 0.210001135
plot(mcp_daytime)
plot(mydf_sp_trans, add=TRUE, col= as.numeric(mydf_sp_trans$daytime)+1)

The daytime activity (or better: hiding and resting) range is ten times smaller than the nighttime activity.

mcp_month <- mcp(mydf_sp_trans[,'month'], percent = 95, unout = "km2") ## MCP (95%) 
mcp_month
## Object of class "SpatialPolygonsDataFrame" (package sp):
## 
## Number of SpatialPolygons:  4
## 
## Variables measured:
##    id      area
## 1   1 0.5789197
## 10 10 0.2092130
## 11 11 0.4173914
## 12 12 0.4405967

Note that there is an activity range calculated for October (month 10), but the animal was only collared 30.10.! This definitely does not make sense.
Now save the MCPs as ESRI shapefiles.

writeOGR(mcp_daytime, dsn = here("output", "data-proc"), 
         "mcp95_daytime_foxQ", "ESRI Shapefile", overwrite = TRUE)

Kernel utility density

Calculate the density kernel:

kud <- kernelUD(mydf_sp_trans[,'daytime'], h = "href") ## calculate kernel with h = "href"
# the output is an object with lists, accessed via: str(kud)
kernel.area(kud, unout = "km2") ## across all levels 5-95%
##        FALSE.       TRUE.
## 20 0.03523210 0.002032032
## 25 0.04844414 0.004064065
## 30 0.06606019 0.006096097
## 35 0.08367624 0.006096097
## 40 0.10129229 0.008128129
## 45 0.11890834 0.010160162
## 50 0.14533242 0.012192194
## 55 0.16735248 0.014224226
## 60 0.19818057 0.016256259
## 65 0.22460465 0.018288291
## 70 0.25543274 0.022352356
## 75 0.29506885 0.024384388
## 80 0.33470497 0.030480485
## 85 0.38755312 0.040640646
## 90 0.45801733 0.058928937
## 95 0.56371363 0.113793810
kud90 <- getverticeshr(kud, percent = 90) ## this creates the spatial object, units m2 
gArea(kud90, byid = TRUE)/ 1e6            ## units km2
##      FALSE       TRUE 
## 0.46582269 0.05652516

Make a plot and save it as spatial polygon:

# plot only the nighttime
image(kud[[1]], col = viridis(100, direction = -1))
xyz <- as.image.SpatialGridDataFrame(kud[[1]])
contour(xyz, add = TRUE)
points(mydf_sp_trans, cex = 0.01, col = "red")

# save it as shapefile - mind - here we can only save one contour line, 
# in this case the 90% kernel:
writeOGR(kud90, dsn = here("output","data-proc"), 
         layer = "kernel_ud90", "ESRI Shapefile", overwrite = TRUE)
# make the crosscheck: load and plot
kernel_ud90_sf <- st_read(dsn = here("output", "data-proc"), layer = "kernel_ud90")
## Reading layer `kernel_ud90' from data source 
##   `C:\Users\DataVizard\Google Drive\Work\Eco\Projects\2022_BioDivCourse\Course4_MoveQ\output\data-proc' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 4556435 ymin: 3270073 xmax: 4557503 ymax: 3271068
## CRS:           unknown
mcp95_sf       <- st_read(dsn = here("output", "data-proc"), layer = "mcp95_daytime_foxQ")
## Reading layer `mcp95_daytime_foxQ' from data source 
##   `C:\Users\DataVizard\Google Drive\Work\Eco\Projects\2022_BioDivCourse\Course4_MoveQ\output\data-proc' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 4556264 ymin: 3270164 xmax: 4557405 ymax: 3271116
## CRS:           unknown
# note: the CRS is lost and not defined! So assign it again:
st_crs(kernel_ud90_sf) <- 3035
st_crs(mcp95_sf)       <- 3035
# TODO eigentlich wollte ich alles übereinander plotten, also mcp, kernel, punkte
## CED: Da implizit im `geom_*` layer aes() als mapping input erwartet wird, man man explizit `data` adressieren, wenn man einen weiteren Datensatz nutzen möchte

ggplot() +
  geom_sf(data = kernel_ud90_sf, aes(fill = id)) +
  geom_sf(data = mydf_sf_trans, size = 1.5, alpha = .1, color = "red3") +
  scale_fill_manual(values = c("grey90", "grey70"), guide = "none") +
  theme_bw(base_size = 16)

Do the following only on a fast computer:

##### der krass coole 3D plot

# install.packages("devtools")
# devtools::install_github("ropensci/plotly")
library(plotly)

my_kud <- kud[[1]] ## or choose kud[[2]]

xy <- coordinates(my_kud)
z <- my_kud@data$ud
df_kud <- data.frame(x = xy[,1], y = xy[,2], z = z)
#persp(x=df_kud$x, y=df_kud$y, z= df_kud$z)
#plot3d(x=xy[,1],y=xy[,2],z=z)

my_kud@grid@cells.dim
## Var2 Var1 
##   60   57
r <- raster(ncols = 60, nrows = 58)
coordinates(df_kud) <- ~x + y
r_kud <- rasterize(df_kud, r, "z", fun = max)

myz <- matrix(z, nrow = my_kud@grid@cells.dim[[2]],
              ncol = my_kud@grid@cells.dim[[1]], byrow = FALSE)
# image(myz)
kd.list <- list(x = xy[,1], y = xy[,2], z = myz)
# with(kd.list, plot_ly(x = x, y = y, z = z, type = "surface"))
plot_ly(z = myz, type = "surface") # click into plot and on icon 'turntable rotation' in plot upper right
#########  ende krass cooler plot ########

Nothing makes sense but in the light of environmental information

-> Exercise 4.2

Exercise 4.2

END

Adding the session info can be very helpful when going back to old scripts or using scripts of others:

Session Info 
sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19043)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252   
## [3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
## [5] LC_TIME=C                      
## system code page: 65001
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] plotly_4.10.0       solartime_0.0.3     viridis_0.6.2      
##  [4] viridisLite_0.4.0   adehabitatHR_0.4.19 adehabitatLT_0.3.25
##  [7] CircStats_0.2-6     boot_1.3-28         MASS_7.3-54        
## [10] adehabitatMA_0.3.14 ade4_1.7-18         deldir_1.0-6       
## [13] move_4.1.6          rgdal_1.5-28        geosphere_1.5-14   
## [16] circular_0.4-93     tmap_3.3-2          ggplot2_3.3.5      
## [19] rgeos_0.5-9         raster_3.5-11       sp_1.4-6           
## [22] sf_1.0-5            lubridate_1.8.0     here_1.0.1         
## 
## loaded via a namespace (and not attached):
##  [1] RColorBrewer_1.1-2      httr_1.4.2              rprojroot_2.0.2        
##  [4] tools_4.1.2             bslib_0.3.1             utf8_1.2.2             
##  [7] R6_2.5.1                KernSmooth_2.23-20      lazyeval_0.2.2         
## [10] DBI_1.1.2               colorspace_2.0-2        withr_2.4.3            
## [13] gridExtra_2.3           tidyselect_1.1.1        leaflet_2.0.4.1        
## [16] compiler_4.1.2          leafem_0.1.6            textshaping_0.3.6      
## [19] xml2_1.3.2              labeling_0.4.2          bookdown_0.24          
## [22] sass_0.4.0              scales_1.1.1            classInt_0.4-3         
## [25] mvtnorm_1.1-3           proxy_0.4-26            systemfonts_1.0.3      
## [28] stringr_1.4.0           digest_0.6.29           rmarkdown_2.11         
## [31] base64enc_0.1-3         dichromat_2.0-0         pkgconfig_2.0.3        
## [34] htmltools_0.5.2         highr_0.9               fastmap_1.1.0          
## [37] htmlwidgets_1.5.4       rlang_0.4.12            farver_2.1.0           
## [40] jquerylib_0.1.4         generics_0.1.1          jsonlite_1.7.2         
## [43] crosstalk_1.2.0         dplyr_1.0.7             magrittr_2.0.1         
## [46] s2_1.0.7                Rcpp_1.0.7              munsell_0.5.0          
## [49] fansi_0.5.0             abind_1.4-5             lifecycle_1.0.1        
## [52] terra_1.4-22            stringi_1.7.5           leafsync_0.1.0         
## [55] yaml_2.2.1              tmaptools_3.1-1         grid_4.1.2             
## [58] parallel_4.1.2          crayon_1.4.2            lattice_0.20-45        
## [61] stars_0.5-5             knitr_1.36              pillar_1.6.4           
## [64] codetools_0.2-18        wk_0.6.0                XML_3.99-0.8           
## [67] glue_1.4.2              evaluate_0.14           leaflet.providers_1.9.0
## [70] data.table_1.14.2       png_0.1-7               vctrs_0.3.8            
## [73] rmdformats_1.0.3        tidyr_1.1.4             gtable_0.3.0           
## [76] purrr_0.3.4             assertthat_0.2.1        cachem_1.0.6           
## [79] xfun_0.27               lwgeom_0.2-8            e1071_1.7-9            
## [82] ragg_1.1.3              class_7.3-19            tibble_3.1.6           
## [85] memoise_2.0.1           units_0.7-2             ellipsis_0.3.2